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‘-  The  effects  of  the  rotation  stress  mechanism  and  the  penetration  of  short¬ 
wave  radiation  below  the  sea-surface  &r^-ex*m4neTd>in  determining  a  one¬ 
dimensional  equilibrium  mixed  layer  depth.  Starting  with  the  Obukhov-scale 
equilibrium  theory  for  the  surface  ocean  boundary  layer,  a  revised 
equilibrium  theory,  which  includes  rotation  stress  and  radiation  effects,  is 
presented.  This  new  theory  is  applied  using  climatological  boundary 
conditions  for  the  tropical  Atlantic^  ah<T the  results  are  compared  with  the 
observed  climatological  mixed  layer  depth. 

In  general,  the  response  of  the  revised  model  is  an  improvement  over  the 
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Obukhov  theory  alone.  Because  the  quality  of  the  results jp#*-  limited  by 
uncertainties  in  the  boundary  conditions,  no  detailed  evaluation  of  the  model 
response  is  justifiable.  However,  it  is  concluded  here  that  the  physical 
mechanisms  of  rotation  stress  and  penetration  of  radiation  are  important  in 
determining  a  steady-state  equilibrium  depth  of  mixing  for  the  tropical 
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I.  INTRODUCTION 


The  study  of  the  oceanic  turbulent  boundary  layer  is  a  relatively  recent 
field  in  Physical  Oceanography.  Since  the  pioneering  work  of  Kraus  and 
Turner  (1967),  a  significant  effort  has  been  made  in  understanding  and 
modelling  the  dynamics  of  turbulent  mixing  in  the  oceans.  However,  the  basic 
contributions  made  by  those  authors  still  provide  the  theoretical  foundation 
for  much  of  the  current  research. 

The  oceanic  boundary  layer  or  mixed  layer  is  considered  in  this  study  to 
be  the  fully  turbulent  region  that  is  bounded  above  by  the  air-ocean 
interface,  and  it  is  where  the  temperature  and  salinity  are  usually  observed 
to  be  fairly  well  mixed.  Below,  the  mixed  layer  is  assumed  to  be  bounded  by  a 
dynamically  stable  thermocline. 

The  study  of  the  top  few  tens  of  meters  of  the  ocean  is  of  considerable 
scientific  interest.  Ekman  pumping  effects,  originating  in  the  mixed  layer  by 
the  action  of  wind  forcing,  considerably  influences  the  dynamics  of  the  lower 
levels  and  the  interior  circulation  in  general.  Also,  the  interaction  between 
the  oceanic  and  atmospheric  boundary  layers  is  an  essential  mechanism 
which  must  be  considered  when  making  medium  and  long-range  weather 
forecasts,  since  a  large  part  of  the  atmospheric  energy  supply  comes  from 
the  heat  exchanged  between  these  layers  .  Additionally, this  is  the  region  of 
primary  biological  productivity,  which  is  of  significant  economic  importance. 
An  important  military  application  is,  on  the  other  hand,  the  modelling  of 
accoustic  propagation  in  the  oceans. 


Host  of  the  physics  behind  existing  models  of  the  oceanic  mixed  layer  are 
based  on  the  flux  form  of  the  Navier-Stokes  equations  of  motion,  with  the 
Bousslnesq  approximation.  One-dimensional  versions,  like  the  one  which  is 
treated  in  this  study,  further  assume  the  horizontal  gradients  of  the  mean 
fields  to  be  negligible.  For  a  steady-state  situation,  this  will  lead  to  the 
basic  energy  balance,  as  stated  by  Kraus  and  Turner  (1967),  between  the 
work  done  by  the  wind  stress  and  the  surface  heating,  in  defining  an 
equilibrium  depth  of  mixing.  However,  a  usually  neglected  term  in  the 
turbulent  kinetic  energy  budgets,  the  planetary  rotation  term,  was  recently 
examined  by  Garwood  et  at.  (1985a).  They  suggest  a  new  formulation  for  the 
equilibrium  mixed  layer  depth,  which  basically  describes  the  interaction 
between  the  zonal  Reynolds  stresses  and  the  northward  component  of  the 
planetary  rotation.  This  new  formulation  might  partially  explain  the  zonal 
variation  of  the  observed  mixed  layer  depth  In  the  Tropical  Pacific  and 
Atlantic  Oceans. 

The  goal  of  this  study  Is  to  apply  the  Kraus  and  Turner  (1967)  and  the 
Garwood  et  ah  (1985b)  one-dimensional,  steady-state  models  on  a  relatively 
broad  grid  with  realistic  boundary  conditions  in  the  Equatorial  and  Tropical 
Atlantic  Ocean,  and  to  compare  the  resulting  diagnostic  mixed  layer  depths 
with  the  observed  mixed  layer  depths.  In  this  way,  insight  may  be  gained  into 
the  relative  importance  of  the  physical  mechanisms  involved.  Also,  we 
should  be  able  to  Identify  the  regions  and  time  of  the  year  where  the  turbulent 
boundary  layer  can  be  represented  realistically  by  such  simple  steady-state 
equilibrium  models. 
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II.  THE  DYNAMICS  OF  MIXING 


A.  THE  TURBULENT  KINETIC  ENERGY  BUDGET 

In  the  presence  of  a  constant  wind  stress  and  a  downward  surface  heat 
flux,  the  steady-state,  one-dimentional  turbulent  kinetic  energy  (TKE) 
budget  is,  assuming  no  vertical  advection  : 

0=d/dt  [E/2 ]= 

-u'w'  du/dz  -v'w'  dv/dz  +  b'w'-d/dz  [w'(E/2)+  w'p'/p]  -  d,  (2.1) 

where  E/2  is  the  TKE,  -u'w'  du/dz  -v’w*  dv/dz  represents  the  production  of 
kinetic  energy  associated  with  the  vertical  shear  in  the  horizontal  flow 
induced  by  the  wind  stress,  b'w'  is  the  buoyancy  flux,  -d/dz  [w(E/2)+w’p7p] 
is  the  turbulent  diffusion  and  d  is  the  dissipation.  With  no  loss  of  generality, 
we  shall  assume  the  x-axis  to  be  oriented  downstream,  sothatv=0. 

1 .  Shear  Production  of  Turbulence 

Let  u(z)  be  some  monotonic  function  of  z  between  the  surface  ,  z=0, 
and  the  bottom  of  the  mixed  layer,  z=-h.  For  the  sake  of  simplicity,  it  is 
further  assumed  that  u(-h)=Q,  since  no  entrainement  is  allowed. 

Consider  now  (Fig.  la)  that  some  particle  initially  at  1  is  displaced  to 
some  position  2  by  means  of  a  positive  turbulent  vertical  velocity  w'. 
Because  du/dz±0,  the  particle  will  acquire  a  negative  turbulent  horizontal 
velocity  u'  and,  then,  u'w'  will  be  less  than  zero  at  2.  Alternatively,  a 
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particle  initially  at  1  f  which  is  displaced  to  3  by  means  of  a  negative  vertical 
turbulent  velocity  w',  will  acquire  a  positive  horizontal  turbulent  velocity  u' 
and,  again,  u‘w'  will  be  negative.  In  both  situations,  the  shear  production 
term  -u'w'  du/dz  will  be  positive,  and  when  this  physical  process  Is 
averaged  over  the  whole  turbulent  layer,  it  will  tend  to  increase  the  TKE  at 
the  expense  of  the  mean  kinetic  energy.  Formally, 

-uV  3u/3z  >  0. 

Evaluated  over  the  mixed  layer,  the  shear  production  will  be*. 

o  ____  _ 

G  =  |  -u'w'  du/dz  dz  =  <-u’w‘>  [u(0)-u(-h)],  (2.2) 

-h 

where  <  >  represents  a  vertical  average.  At  the  bottom  of  the  mixed  layer, 
and  under  the  assumption  of  no  entrainement,  -uVl-h)  is  zero.  Then, 
disregarding  nondlmenslonal  constants,  <-t7w'>  Is  of  order  -iTw'(O).  Also, 
we  have  already  assumed  u(-h)  to  be  zero.  Thus, 

G  =  -CTw*(0)  u(0).  (2.3) 

Furthermore,  the  Reynold  stress  -uV(0)  can  be  written 


where  p  is  a  typical  sea  water  density  and  t  is  the  surface  wind  stress. 
Finally,  (2.3)  is  rewritten  as 


G  =  u(0)|t|/p  . 


(2.5) 


2.  The  Buovancv  Flux 

Under  the  influence  of  surface  heating,  there  will  be  a  tendency  for 
stratification,  and  mixing  will  only  occur  in  the  presence  of  wind  stirring.  In 
this  case,  the  less  dense  particles  near  the  surface  will  mix  with  the 


Fig.  1.  (a)  Shear  production  of  turbulent  kinetic  energy,  in  the  presence  of 
wind  stress,  (b)  Buoyant  damping  of  turbulent  kinetic  energy,  in  the 
presence  of  surface  heating. 
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underlying  water  only  at  the  expense  of  the  wind-generated  turbulence,  thus 
converting  TKE  into  potential  energy  (Fig.  lb).  This  means  that,  on  the 
average,  the  buoyancy  flux  -bV'  will  be  a  sink  term  in  the  TKE  budget.  The 
equilibrium  between  the  effects  of  wind  mixing  and  surface  heating  is  indeed 
the  fundamental  energy  balance  for  any  steady-state,  one-dlmentlonal  mixed 

layer  model,  as  long  as  advectlon  Is  ntjected. 

To  relate  the  buoyancy  flux  to  the  surface  heat  flux  Q0  we  apply  the 

linearized  equation  of  state  and  obtain: 

b'w'sagTV,  (2.6) 

where  a  is  the  thermal  expansion  coefficient,  g  is  the  acceleration  of 
gravity  and  TV  is  the  temperature  flux.  It  can  further  be  shown  that, 
assuming  vertical  homogeneity  In  the  mean  temperature  T,  TV  will  be  linear 
over  the  mixed  layer,  and 


TV'(O)  »  -Q0/(  p  Cp),  (2.7) 

where  Q0,  the  surface  heat  flux  is  positive  downward. 

Integrating  the  buoyancy  flux  over  the  mixed  layer, 

o o  _  _____  _____ 

B2  J  bV  dz  =  ag  J  TV  dr  =  agh  [TV’(0)+TV(-h)]/2,  (2.8) 

-h  -h 

where  h  Is  the  mixed  layer  depth.  Because  no  entrainment  is  assumed, 
fw'(-h)=0.  Thus, 
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B  -  -  agh  Q0/(2pCp). 


(2.9) 


This  simple  model  for  the  vertically  integrated  buoyancy  flux  also 
requires  the  assumption  that  the  net  solar  radiation  be  completly  absorbed  at 
the  surface.  We  will  see  in  the  following  sections  that  this  assumption  is 
incorrect  when  the  vertical  scale  of  penetration  of  radiation  below  the 
surface  is  comparable  to  the  vertical  scale  of  turbulent  mixing.  The  effect  of 
the  penetration  of  radiation  is  expected  to  be  relevant  in  optically  clear 
oceanic  waters  and  in  summertime  conditions,  when  the  mixed  layer  is 
relatively  shallow.  In  those  situations,  the  net  downward  buoyancy  flux  may 
be  significantly  reduced,  causing  the  resultant  depth  of  mixing  to  be 
significantly  increased. 

3.  Turbulent  Diffusion 

Although  the  diffusion  of  turbulence  Is  an  essential  mechanism  in 
transporting  TKE  downwards  from  the  upper  levels,  It  does  not  constitute,  by 
itself,  a  source  or  sink  of  turbulence,  which  means  that  its  net  value  over 
the  mixed  layer  must  vanish.  Formally, 

I  -d/dz  lw'(E/2)+w’p7p]  dz  =  0.  (2.10) 

-h 

4.  PlSStBLflllQn 

Early  mixed  layer  models,  including  the  Kraus  and  Turner  (1967) 
model,  recognized  the  existence  of  dissipation.  However,  they  considered  its 
influence  to  be  negligible,  given  the  relatively  large  length  scale  of  turbulent 


mixing.  Indeed,  a  less  careful  scale  analysis  of  the  TKE  budget  seems  to 
support  the  approximation.  More  recent  models,  however,  include 
dissipation  In  the  mechanical  energy  budget,  either  as  a  fraction  of  the  wlnd- 
-stress  production  or  as  a  fraction  of  the  net  TKE  production.  Although 
viscous  dissipation  may  be  neglected  when  dealing  with  turbulent  eddies  of 
length  scales  of  order  of  the  mixed  layer  depth,  the  dynamic  interaction 
between  those  eddies  produces  other  eddies  of  smaller  and  smaller 
dimensions,  where  the  role  of  dissipation  may  finally  dominate.  This  fact  is 
easily  understood  with  the  analysis  of  a  spectrum  of  turbulence,  where  it  can 
be  verified  that  the  mechanical  energy  drops  rapidly  above  a  limiting 
frequency  (Tennekes  and  Lumley,  1962). 

If  such  a  dissipation  mechanism  Is  Included  in  the  TKE  budget  by 
parameterization,  for  example  as  a  fraction  of  the  TKE  production,  we  will 
see  in  the  following  sections  that  the  Kraus  and  Turner  model  is  still  usable, 
since  no  aditlonal  degree  of  freedom  Is  required.  That  was  the  approach  of 
Gelsler  and  Kraus  (1969),  Nlller  (1975)  and  others. 

For  now,  we  will  just  define  the  vertically  Integrated  dissipation  D  as: 

0  »  f  d  dz .  (2.11) 

-h 


B.  EQUILIBRIUM  DEPTH  OF  MIXING 

From  the  results  of  the  previous  section,  the  vertically  integrated  TKE 
budget  can  be  written  as 
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where 

and 


B  +  G  -  D  =  0, 


(2.12) 


B  -  -  agh  G0/(pCp) 
G  =  u(0)W/p  . 
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Fig.  2.  The  one-dimensional  mixed  layer  equilibrium.  The  'redistribution* 
effects  include  the  pressure  redistribution  and  the  rotation  stress 
mechanism. 
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Equation  (2.12)  illustrates  the  basic  steady-state  balance  between  the  work 
done  by  the  wind  stress  and  the  buoyant  damping  (Fig.  2).  Note  that,  because 
of  turbulent  mixing,  the  mean  flow  will  tend  to  be  vertically  homogeneous, 
except  near  the  surface,  where  turbulence  is  generated,  and  near  the  bottom 
of  the  mixed  layer,  where  It  falls  to  zero.  Also,  the  mean  temperature 
profile  will  tend  to  be  homogeneous.  Thus,  the  mixed  layer  potential  energy 
will  approach  Its  theoretical  maximum  value. 

In  Fig.  2  an  adltlonal  mechanism  is  Illustrated.  This  processs,  called 
‘redistribution’,  does  not  appear  in  the  TKE  budget,  since  its  constituent 
terms  sum  to  zero  when  the  three  spatial  components  of  the  kinetic  energy 
budget  are  added.  Redistribution  consist  of  two  parts:  the  pressure 
redistribution  and  the  planetary  rotation.  Their  roles  are  to  exchange  TKE 
among  the  x-y-z  directions,  tending  to  restore  the  system's  isotropy.  In 
section  III,  the  planetary  rotation  term  will  be  shown  to  have  a  fundamental 
role  in  determining  the  equilibrium  mixed  layer  depth  for  the  Garwood  et  a/. 
(1985a)  steady-state  model, 
i.  The  Obukhov  Length  Scale 

We  have  seen  that,  if  dissipation  Is  parameterized  in  terms  of  the 
other  basic  turbulent  variables,  the  steadiness  of  the  mixed  layer  depth  will 
depend  on  the  surface  boundary  conditions,  l.e.,  the  surface  wind  stress  and 
the  surface  heat  flux.  This  Is  equivalent  to  setting  G/B  =  constant. 

We  will  now  define  a  velocity  scale  u*  ,  which  is  representative  of  the 

typical  turbulent  speed  of  the  largest  eddies  in  a  wind-driven  mixed  layer: 

u*2-|T|/p.  (2.13) 
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Because  u(0)  depends  linearly  on  the  wind  speed  and,  thus,  on  the  surface 
wind  stress,  the  shear  production  term  G  In  equation  (2.5)  can  be  scaled  as 

G«-u*3.  (2.14) 


Combining  (2.9)  and  (2.12-14),  and  setting  -G/B  =1,  the  Obukhov 
length  scale  is  defined  as 


L  =  h  =  2  p  Cp  u*3  /  (a  g  Q0).  (2.15) 

Physically,  L  represents  the  scale  of  maximum  depth  of  turbulent  mixing, 
given  the  surface  boundary  conditions  of  wind  stress  and  heat  flux,  for  an 
equilibrium  mixed  layer  with  no  advectlon  present. 

2.  im  Kraus  and  Turner  Steady-State  Model 

Based  on  equation  (2.12),  Kraus  and  Turner  (1967)  considered  the 
dissipation  to  be  negligible,  so  that 

B  +  G  =  0.  (2.16) 


Also,  the  vertically  Integrated  shear  production  G  can  be  written 


or 


G  »  u*3  ■  (|it|/p)3/2 
G-  (pa  CDu2/p)3/2  , 


(2.17) 

(2.18) 


where  Cp  is  the  drag  coefficient,  pa  is  the  air  density  and  U  is  the  wind 

speed  measured  at  10  meters  above  the  surface.  Substituting  (2.9)  and 
(2.18)  into  (2.16),  the  equilibrium  depth  of  mixing  is  written 

h  =  Z)  (  p  Cp/  a  g  Q0  )(  Pa  CD  u2/p)3/2  ,  (2. 19) 

where  Cj  is  a  dimensionless  constant  of  order  1 . 

More  recent  versions  of  the  Kraus  and  Turner  model  assume  the 
dissipation  D  to  be  proportional  to  the  total  TKE  production  B+G,  so  that 
equation  (2.12)  is  written  as 

(K+1)(B  +  G)  =  0,  (2.20) 

where  K  is  a  dimensionless  constant.  As  stated  earlier,  this  will  lead  to  the 
same  basic  result  of  equation  (2. 19),  since  no  aditiona)  degree  of  freedom  is 
included. 


III.  THE  REVISED  STEAD'  STATE  MODEL 


As  noted  in  section  I,  the  dynamics  of  the  oceanic  mixed  layer  in  the 
equatorial  and  tropical  regions  are  not  fully  understood.  One  of  the  most 
apparent  but  least  understood  features  is  the  zonal  variation  of  the  mixed 
layer  depth. 

Some  of  the  discrepancies  between  the  observed  mixed  layer  depths  and 
the  predictions  made  by  existing  steady-state  and  one-dimensional  models 
might  be  explained  by  the  exceptionally  strong  vertical  and  horizontal 
advection  which  are  characteristic  of  the  tropical  oceans.  Also,  the  effect  of 
penetration  of  short  wave  radiation  below  the  surface  is  believed  to  have 
some  importance  for  relatively  shallow  mixed  layers,  especially  in  the 
presence  of  large  net  solar  irradiance.  Another  physical  mechanism, 
examined  by  Garwood  et  a/.  (1985a),  involves  the  interaction  between  the 
zonal  surface  wind  stress  and  the  planetary  rotation,  and  it  may  partially 
explain  the  zonal  dependence  of  the  steady-state  mixed  layer  depth. 

It  is  beyond  the  scope  of  this  work  to  formally  treat  the  effects  of 
advection  on  the  dynamics  of  the  mixed  layer,  and  we  shall  restrict  our  study 
to  steady-state  situations.  However, the  other  two  mechanisms  referred  to 
(penetration  of  radiation  and  planetary  rotation)  will  be  included  in  the  model 
derived  in  the  following  sections. 
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A.  EFFECT  OF  THE  PENETRATING  RADIATION  ON  THE  BUOYANCY  FLUX 

The  simplest  model  for  the  penetration  of  radiation  below  the  sea-surface 
assumes  that  the  net  long-wave  solar  irradlance  (  Infra-red  and  red  )  is 
absorbed  at  the  surface,  and  the  short-wave  irradlance  decays  exponentially 
with  depth.  Recent  models  (  Paulson  and  Simpson,  1977  ;  Simpson  and  Dickey, 
1980)  also  allow  for  the  penetration  of  long-wave  radiation.  In  both 
approaches,  the  predicted  vertical  profile  of  absorbed  irradlance  is  sensitive 
to  the  values  of  two  empirical  parameters:  the  vertical  scale  of  penetration, 
which  mainly  depends  on  the  turbidity  of  the  water,  and  the  short-wave 
(blue-green)  fraction  of  the  net  solar  irradlance  incident  on  the  surface. 
Because  there  is  only  a  limited  amount  of  data  on  the  optical  properties  of  the 
oceans,  the  values  of  these  empirical  parameters  will  either  have  to  be 
assumed  or  adjusted  as  tunning  constants.  For  that  reason,  and  also  because 
a  steady-state  model  will  be  applied  to  large  oceanic  regions  having  variable 
optical  properties,  the  simplest  model  will  be  used,  i.,e.,  only  the  short¬ 
wave  radiation  will  be  assumed  to  penetrate  below  the  surface.  The  vertical 
profile  of  the  penetrating  irradlance  will  have  therefore  the  following 
exponential  form 


Q2(z)  -  RQS  ez/X  •  (3.D 

Here,  the  downward  short-wave  irradlance,  Q2U),  is  the  radiant  flux  density 

on  a  horizontal  surface  due  to  contributions  from  the  entire  upward 
hemisphere;  Qs  is  the  net  solar  irradiance  at  the  surface,  R  is  the  short-wave 

fraction  of  the  net  solar  Irradlance  (  blue-green  ),  z  Is  the  vertical 
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coordinate,  positive  upward;  and  X  is  the  penetration  scale,  assumed 
constant  with  depth. 

Alternatively,  the  long-wave  solar  irradtance,  assumed  to  be  totally 
absorbed  at  the  surface,  is  given  by 


Q|  *  Qg  —  R  Qg  ,  (3.Z) 

where  Q0  *  Qs-  Q&-  Q^-  Qe  is  the  net  solar  irradiance  at  the  surface,  minus 

the  long-wave  back  radiation,  minus  the  sensible  heat  flux,  minus  the  latent 
heat  flux. 

Assuming  horizontal  homogeneity  of  the  variables,  it  can  be  shown  that 
the  effective  heat  flux  Q  contributing  to  the  temperature  flux  over  the  mixed 
layer  will  have  the  form  (  Garwood,  1977  ) 

Q  -  Qp  f°l  dQ2/  dz  -  2/h  fdQe/ds  ds  ]  dz  (3.3) 

or  Q*Q0-RQsl2X/h-  e  _h/X  (  1  +  2X/h  )).  (3.4) 

For  deep  mixed  layers  (  h/X  »1)  or,  puting  R  =  0,  expression  (3.4) 
reduces  to  Q  »  Q0,  which  corresponds  to  the  net  heat  flux  being  absorbed  at 

the  surface.  For  shallow  mixed  layers  (  h/X«1),  expression  (3.4)  reduces 
to  Q  *  Q0  -  R  Qs.  Under  these  circumstancies  the  effect  of  the  penetrating 

radiation  on  the  effective  heat  flux  is  maximum. 

It  is  also  apparent  from  expression  (3.4)  that,  In  summertime  conditions 
(  small  h,  large  Qs  )  and  in  clear  oceanic  waters  (  large  X),  Q  can  be 
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significantly  smaller  than  Q0,  thus  reducing  the  buoyancy  flux  and  increasing 

the  vertical  extent  of  the  wind  generated  turbulence. 

From  expression  (2.9),  the  effective  vertically  integrated  buoyancy  flux 
will  then  be 


B  =  -  ag  h  Q  /(  2pCp  ), 


(3.5) 


where  Q  is  given  by  (3.4). 


A.  THE  ROTATION  STRESS  MECHANISM 

The  steady-state  turbulent  kinetic  energy  budgets  for  the  three  spatial 
components  (  x  eastward,  y  northward,  z  upward  ),  vertically  integrated 
over  the  mixed  layer  are: 

0  ■  Gx  +  Px  +  Qz  h  UV(0)  +  0 y  hTx/p  -  D/3  (3.6a) 

0  =  Gy  +  Py-GzhuV(0)  -D/3  (3.6b) 

0  =  B  +  Pz  -  Oy  hTx/p  —  D/3.  (3.6c) 

Here  Gx  and  Gy  are  the  horizontal  components  of  the  shear  production,  B  Is 
the  vertically  Integrated  buoyancy  flux,  the  P\'s  are  the  components  of  the 
vertically  integrated  turbulent  transport,  and  tx  Is  the  zonal  component  of 
the  surface  wind-stress.  The  flj's  are  the  two  spatial  components  of  the 
planetary  rotation  and  D  is  the  dissipation,  assumed  to  be  equally  partitioned 
among  components  (local  isotropy  assumed).  Since  the  transport  and 
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planetary  rotation  terms  vanish,  the  sum  of  the  three  above  equations  simply 
leads  to : 

B  +  G-D  =  0.  (2.12) 

An  aditlonal  process,  not  apparent  in  (2.12),  Is  evident  in  equations 

(3.6).  This  is  the  exchange  between  horizontal  and  vertical  turbulent  kinetic 
energy,  in  the  presence  of  the  zonal  wind  stress  tx  and  the  meridional 

component  Qy  of  the  planetary  rotation.  Note  that  the  sign  of  tx  will 
determine  whether  the  term  Qyhtx/p  is  a  source  or  sink  of  horizontal 
turbulent  kinetic  energy.  With  easterly  winds  (  tx  <  0  ),  the  horizontal 
kinetic  energy  is  expected  to  be  enhanced  at  the  expense  of  the  vertical 
kinetic  energy,  with  a  consequent  shallowing  of  the  mixed  layer.  With 
westerly  winds  (  tx  >  0  ),  the  exchange  between  the  vertical  and  horizontal 

components  will  have  the  opposite  sign,  leading  to  an  increase  of  the  vertical 
extent  of  turbulent  mixing.  This  mechanism,  usually  overlooked  in  planetary 
boundary  layer  models,  was  first  examined  by  Garwood  et  a/.  (1985b),  in  an 
attempt  to  explain  the  unusually  deep  extent  of  mixing  in  the  central  and 
western  Pacific. 

If  the  dissipation  Is  parameterized  In  terms  of  the  vertical  average  of  the 
total  turbulent  kinetic  energy,  1 . e, 

0  =  <  E  >3^2 f  (3.7) 
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where  rri|  and  subsequent  mfs  are  dimensionless  constants,  (2.12)  can  be 


written 


G  +  B- mi  <E  >3/2  =  0  .  (3.8) 

Also,  the  vertical  component  Pz  of  the  turbulent  transport  can  be 

parameterized  in  terms  of  the  vertical  average  of  the  TKE.  According  to  Rotta 
(1951 ),  a  first  order  approximation  to  the  transport  term  Pz  is 

Pz  =  m2<E>3/2.  (3.9) 

Substituting  (3.9)  into  (3.6c)  and  combining  with  (3.8),  will  yield  a 
system  of  two  equations  which  can  be  solved  for  the  mixed  layer  depth  h  and 
the  vertical  average  of  the  turbulent  kinetic  energy  <  E  >,  given  the  surface 
boundary  conditions  and  the  meridional  component  of  the  planetary  rotation. 
The  predicted  value  of  h  will  then  reflect  the  interaction  between  the  zonal 
wind  stress  and  the  Earth  rotation. 

C.  EQUILIBRIUM  DEPTH  OF  MIXING 
We  define 


L  =  2pCp  u.3/(  ag  Q  ) 

(3.10) 

G  =  m3  u*3 

(3.11) 

and 

$  =  fly"txE-p^(  dg  Q), 

(3.12) 
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where  the  dimensionless  parameter  ♦  is  the  quotient  between  the  conversion 
of  vertical  turbulent  kinetic  energy  by  means  of  the  rotation  mechanism  and 
the  buoyant  damping  of  vertical  turbulent  kinetic  energy.  Then,  a  diagnostic 
equation  for  the  mixed  layer  depth  can  be  derived: 

h  *  C1  L  /(!+  C2  ♦),  (3.13) 

where  the  non-dimensional  constants  and  C2  are  related  to  the  m{s: 

Cj  =  m3(m2/mpl/3)/(m2/mj+2/3)  (3.M) 

and  C2  -  2/(m2/m1+2/3).  (3.15) 

Note  that  the  effective  surface  heat  flux  Q  In  equations  (3. 10)  and  (3. 12) 
reflects  the  penetration  of  radiation  below  the  surface,  as  modelled  by 
equation  (3.4).  Thus  the  diagnostic  value  of  h  will  be  dependent  on  the  values 
of  four  tunning  parameters:  the  Cj's,  the  factor  R  and  the  vertical  scale  of 

penetration  X.  Because  Q  is  a  relatively  complex  function  of  h,  a  solution  for 
h  from  (3.4)  and  (3.13)  is  only  possible  by  iterative  techniques.  An 
alternative  method,  which  significantly  reduces  the  computational  effort,  is 
to  substitute  the  Obukhov  length  scale  L  for  h  In  equation  (3.4),  and  then 
solve  for  Q.  This  simplification  is  justified  and  will  be  utilized  in  the  next 
sections,  since  the  primary  goal  of  this  study  is  to  examine  the  gross 
characteristics  of  the  solution  space. 

Substituting  (3.13)  back  into  (3.7),  the  following  expression  for  the 
vertically  integrated  dissipation  is  obtained: 
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D«  m,  <  E  >5/2  -  C2  m3  u»*  (l+ZD/UU+C?*)]. 


(3.16) 


Physically,  the  dissipation  must  always  be  non-negative,  which  yields 


1+2*  iO 


♦  i  -1/2. 


(3.17) 


Fig.  3.  The  dimensionless  depth  scale  h*  plotted  against  the  rotation  stress 
variable  ®. 
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For  ♦  <  -1/2,  and  In  absence  of  advection,  no  steady-state  situation  is 
possible  and  the  mixed  layer  will  continually  entrain.  Also,  noting  that,  in 

equation  (3.13),  the  denominator  must  be  positive,  the  limiting  case 
expressed  by  (3.17)  suggests  that  C2  <  2  for  all  situations  (since  C2  is 

assumed  to  be  a  constant).  Zeman  and  Tennekes  (1975)  obtained 
m2/mj~1/2,  which  gives  C2 ~  12/7.  Using  the  suggested  value,  a 

nondimensional  mixed  layer  depth  is  given  by  (  Fig.  3  ) 

h*  *  h/  CtL  *  1/(  1  +  12#/  7  ).  (3.18) 

Note  that  to  the  limiting  value  ♦  =  -1/2  corresponds  a  theoretical  maximum  of 
h*  =  7. 
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IV.  THE  EQUILIBRIUM  MIXED  LAYER  DEPTH  IN  THE  TROPICAL  ATLANTIC 


A.  THE  BOUNDARY  CONDITIONS 

In  the  following  sections,  the  method  of  application  of  the  steady-state 
mixed  layer  model  to  an  oceanic  macro-scale  domain  Is  presented.  The  goal 
of  this  study  is  to  gain  insight  into  the  relative  importance  of  the  physical 
mechanisms  Involved,  as  well  as  to  identify  the  regions  and  the  time  of  the 
year  where  the  observed  turbulent  layer  can  be  represented  realistically  by 
such  a  simple  steady-state  model.  This  goal  is  limited  to  some  degree  by  the 
availability  and  suitability  of  experimental  observations  of  the  boundary 
conditions.  Although  a  significant  effort  has  been  made  in  the  last  decade  to 
obtain  accurate  and  detailed  measurements  of  the  surface  boundary 
conditions,  such  data  are  frequently  confined  to  limited  oceanic  areas.  The 
available  climatological  data  for  large  oceanic  areas  are  often  unreliable, 
sparse  and  distributed  Irregularly,  depending  on  the  usual  routes  followed  by 
mariners.  Given  these  limitations,  the  tropical  Atlantic  Ocean  is  the  best 
covered  ocean  basin. 

Except  for  the  observed  surface  wind-stress,  where  the  digitized  data  of 
Hellerman  and  Rosensteln  (1983)  were  used,  all  the  other  data  were 
Interpolated  from  climatological  atlases.  The  area  chosen  was  the  tropical 
and  equatorial  Atlantic  between  6  degrees  south  and  30  degrees  north. 

1.  The  Observed  nixed  Laver  Depth 

The  atlas  of  Robinson  et  at.  (1979)  was  used  to  obtain  the  observed 
monthly  values  of  the  mean  depths  to  the  top  of  the  thermocline  for  the 
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equatorial  and  tropical  Atlantic.  The  top  of  the  thermocline  is  therein  defined 
as  the  depth  at  which  the  temperature  is  1 . 1  degrees  centigrade  less  than  he 
surface  temperature.  Although  this  definition  differs  from  the  usual 
definition  of  'mixed  layer  depth’  (the  vertical  extent  of  turbulent  mixing),  the 
difference  should  not  be  significant  for  the  present  study,  provided  density  is 
mostly  dependent  upon  temperature,  and  the  boundary  layer  Is  well  mixed. 
For  the  cases  where  the  mixed  layer  mean  temperature  decreases 
significantly  with  depth,  the  above  criterion  might  underestimate  the  vertical 
extent  of  mixing. 


Fig.  4.  Zonal  variation  of  the  mean  mixed  layer  depth  for  the  month  of 
September  (equator).  The  curve  is  the  result  of  an  adjusted  cubic  spline  fit 
over  the  observed  values,  represented  in  the  figure  by  circles. 


Because  the  contours  displayed  In  the  atlas  have  a  relatively  large 
Increment  of  15  meters,  a  cubic  spline  fit  was  made  for  each  parallel  of 
latitude,  and  the  resulting  curve  was  further  adjusted  so  that  the  interpolated 
values  for  each  interval  would  lie  between  extremes.  The  resulting  maximum 
graphical  error  Is  estimated  to  be  ±10  meters.  No  indication  fs  given  by  the 
authors  about  the  estimated  errors  for  their  analtzed  depths. 

Fig.  4  shows  the  zonal  equatorial  profile  of  the  mean  mixed  layer 
depth  for  the  month  of  September,  as  computed  by  the  above  technique. 

2.  The  Surface  Heat  Flux 

The  heat  budget  atlas  of  Hastenrath  and  Lamb  (1978)  was  used  to 
provide  best  estimates  of  the  net  surface  heat  flux  and  the  net  surface  solar 
irradlance  for  the  area.  For  the  assumptions  made  and  the  calculation 
details,  the  reader  is  referred  to  the  introductory  part  of  the  atlas. 

The  contour  interval  for  the  surface  heat  flux  is  40  W/m2,  which 
constitutes  a  low  resolution.  In  order  to  minimize  the  resulting  errors,  a 
cubic  spline  surface  was  first  fitted  over  a  grid  of  digitized  values,  and  the 
resulting  interpolated  values  were  further  adjusted  for  each  parallel  of 
latitude.  Thus,  for  each  interpolation  interval,  the  calculated  values  would 
lie  between  extremes.  The  considerable  effort  required  for  this  process  is 

well  Justified  because  of  the  sensitivity  of  the  diagnostic  equations  to  small 
deviations  in  the  values  of  Q0,  specially  when  Q0  is  small. 

The  relative  and  absolute  errors  in  the  net  surface  heat  flux,  as 
indicated  by  the  authors,  are  estimated  at  less  than  10  and  more  than 
20  W/m2,  respectively.  The  absolute  error  in  interpolated  values  is 
estimated  to  be  less  than  20  W/m2. 
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As  an  example,  Fig.  7a  shows  the  zonal  equatorial  profile  of  the  mean 
surface  heat  flux  for  the  month  of  September,  as  determined  by  the  above 


technique. 

3.  The  Surface  Wind  Stress 

The  surface  wind  stress  data  are  those  of  Hellerman  and  Rosensteln 
(1983),  for  a  mesh  grid  of  2  degrees  of  latitude  by  2  degrees  of  longitude, 
which  is  suitable  for  our  purpose.  The  two  components  are  calculated  by 
means  of  the  bulk  aerodynamic  formulation: 

TX  *  PacD  Ux  (Ux2  +  Uy2  )1/2  (4. 1 ) 

and  Ty  -  paCo  Uy  (Ux2  +  Uy2  )1/2,  (4.2) 

where  Ux  and  Uy  are  the  two  components  of  the  wind  speed  at  10  meters 
above  the  sea-surface.  The  drag  coefficient,  Cp,  depends  upon  wind  speed 

and  stability,  according  to  the  formulation  of  Bunker  (1976).  The  maximum 
standard  error  in  the  surface  wind-stress  values  for  the  area  considered  is 
estimated  to  be  less  than  0.01  N/m2. 

The  wind  stress  data  are  given  for  each  odd  whole  degree  of  latitude 
and  longitude,  and  In  the  present  scheme  we  will  be  working  with  even 
parallels  of  latitude.  Thus  the  original  data  had  to  be  interpolated.  The 
resulting  interpolation  errors  are  considered  to  be  negligible  for  the  purpose 
of  the  present  study. 

4.  Steadv-State  versus  Time-Dependent  Solutions 


Formally,  a  steady-state  mixed  layer  model  should  only  be  applied  to 
a  real  oceanic  situation  when  entrainment  has  stopped  and  the  top  of  the 


thermocllne  has  retreated  to  an  equilibrium  level,  dependent  on  the  balance 
between  the  downward  surface  heat  flux  and  the  surface  wind-stress.  On  a 
seasonal  time  scale,  i.e.,  neglecting  daily  variations,  and  for  mid- 
latitudes,  the  above  situation  is  expected  to  occur  in  spring  and  summer, 
when  the  mean  downward  surface  heat  flux  is  rising  and  the  wind  speed  is 
decreasing.  For  the  equatorial  and  tropical  regions,  however,  this  process 
is  complicated,  since  the  mean  surface  heat  flux  follows  a  complex  annual 
variation,  while  the  wind  speed  has  a  strong  annual  component. 

For  simplicity,  the  diagnostic  mixed  layer  depth  should  be  calculated 
for  the  same  month  throughout  the  chosen  area,  though  the  requirement  of 
steadiness  will  not  be  satisfied  in  some  regions.  A  logical  and  simple  way  to 
solve  this  problem  is  to  plot  the  annual  time  series  of  the  zonal  profiles  of 
the  observed  mixed  layer  depths  for  some  parallels  of  latitude  and  choose 
the  time  of  the  year  which  best  conforms  to  the  condition  of  steadiness  for  all 
latitudes.  Figure  5  shows  these  plots,  as  compiled  and  further  interpolated 
from  Robinsons's  atlas.  The  month  of  September  was  chosen  as  conforming 
best  to  these  conditions.  Except  at  the  equator,  where  the  top  of  the 
thermocline  seems  to  be  still  deepening,  the  other  latitudes  show 
approximate  steady-state  situations,  during  September. 

The  analysis  of  Fig.  5a  reveals  a  strong  annual  component  of  the 
equatorial  mixed  layer  depth,  which  suggests  the  domination  of  the  wind 
stress  in  determining  the  vertical  extent  of  mixing.  In  Fig.  6,  the  annual  time 
series  of  the  net  surface  heat  flux  and  the  wind  speed  are  shown.  Indeed,  a 
qualitative  compalrason  of  the  contour  lines  of  both  figures  (5a  and  6) 
reveals  a  much  stronger  correlation  between  the  mixed  layer  depth  and  the 
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wind  speed  than  between  the  mixed  layer  depth  and  the  surface  heat  flux.  In 
particular,  the  wind  speed  maximum  (August  to  October)  and  minimum 
(March  to  April),  coincide  well  with  the  corresponding  relative  extremes  in 
the  mixed  layer  depth.  Similar  agreement  is  not  apparent  between  the  mixed 
layer  depth  and  the  surface  heat  flux,  as  so  far  as  an  equilibrium  state  is 
concerned. 

On  the  other  hand,  one  might  expect  the  annual  variation  of  the  surface 
heat  flux  to  show  a  marked  response  due  to  the  annual  variation  of  the  Sun’s 
declination,  i.e.,  to  show  a  strong  biannual  component.  However,  the 
exchange  of  heat  between  the  ocean  and  the  atmosphere  is,  in  the  equatorial 
and  tropical  regions,  significantly  Influenced  by  other  variables.  These 
variables,  such  as  cloudiness  and  precipitation,  have  seasonal  variations 
which  are  more  complex  and  difficult  to  explain. 

During  the  months  of  May  to  September, another  Interesting  feature  in 
the  eastern  part  of  the  Equator  is  the  mixed  layer  depth  minimum.  This 
occurs  when  the  surface  heat  flux  is  negative  (upward).  Without  advection, 
an  equilibrium  situation  is  obviously  precluded.  Interestingly,  the  atlas  of 
Robinson  et  a!.  (1979)  shows  that  this  is  also  a  zone  of  strong  equatorial 
upwelling  and  surface  divergence.  Thus,  an  advective  mechanism  may 
largely  Influence  the  dynamics  of  turbulent  mixing  for  that  region. 

From  the  above  discussion  and  also  from  the  knowledge  that  the 
equatorial  and  tropical  regions  are  zones  of  important  oceanic  heat 
divergence,  it  might  be  concluded  that  horizontal  and  vertical  advection  play 
an  important  role  in  determining  the  vertical  extent  of  mixing  at  the  equator. 
As  lateral  transport  in  the  ocean  diverts  much  of  the  heat  exchanged  with  the 
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atmosphere,  this  effect  would  tend  to  spactally  homogenize  the  surface  heat 
flux  available  for  the  vertical  turbulent  processes.  Thus,  the  seasonal 
variations  of  the  mean  mixed  layer  depth  would  primarily  reflect  the 
seasonal  variations  of  the  wind  speed. 

B.  METHOD  OF  MODEL  TUNING 

The  models  described  in  sections  II  and  III  shall  be  applied  to  a  set  of 
some  parallels  of  latitude  In  the  tropical  and  equatorial  Atlantic,  from  6 
degrees  south  to  30  degrees  north.  For  each  parallel,  the  spacing  of 
calculated  values  will  be  1  degree  of  longitude.  The  purpose  of  the  following 
procedure  is  to  identify  the  zones  where  the  observed  mixed  layer  is  well 
represented  by  the  above  steady-state  models,  as  well  as  to  analyse  the 
relative  Influence  of  the  physical  mechanisms  discussed  In  section  III:  the 
penetration  of  radiation  below  the  surface  and  the  rotation  stress 
mechanism.  For  each  parallel  of  latitude,  a  zonal  profile  of  the  following 
nondimensional  depth  will  be  first  calculated: 

H*  =  H  (1+C2  0)/L  ,  (4.3) 

where  H  is  the  observed  mixed  layer  depth  and  C2  =  12/7.  To  calculate  the 
Obukhov  length  scale  L  and  the  rotation  stress  variable  0,  the  penetration 
parameters  will  be  set  to  ,  respectively,  R=0.2  and  A.  =  6  meters.  Note 
that,  for  each  longitude  point,  H*  is  numerically  equal  to  the  constant  C^.  In 

the  ideal  case  of  a  perfect  agreement  between  the  model's  results  and  the 
observations,  the  resulting  curve  of  H*  for  each  latitude  would  be  a  constant 
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of  order  l.  The  departure  of  the  calculated  values  from  that  ideal  solution 
will  constitute  a  method  of  identifying  the  zones  where  steadiness  is  not 
verified  or  advection  plays  an  important  role.  For  the  situations  where  the 
surface  heat  flux  Is  upward,  i.e. ,  where  a  one-dimensional  equilibrium  is  not 
possible,  H*  will  be  set  to  0.  The  described  procedure  also  makes  possible 

the  selection  of  the  areas  where  the  models  will  be  applied.  For  these  areas, 
the  below  defined  diagnostic  mixed  layer  depths  will  then  be  calculated: 

(i)  h0  :  the  equilibrium  mixed  layer  depth,  assuming  R=0,  0^=0.  This 
corresponds  to  the  steady-state  Kraus  and  Turner  diagnostic  depth,  where 

radiation  is  assumed  to  be  totally  absorbed  at  the  surface. 

(ii)  h^  :  the  equilibrium  mixed  layer  depth,  assuming  R=0.2,  X  =  6m, 
C2=0.  This  is  the  same  as  h0,  except  that  radiation  is  assumed  to  penetrate 

below  the  surface,  according  to  the  model  given  by  equation  (3.4). 

(iii)  t>2  *•  the  equilibrium  mixed  layer  depth,  assuming  R=0.2,  X  *  6m  and 
C2*0.  Here,  the  rotation  stress  mechanism  is  also  included,  according  to  the 
formulation  of  Garwood  et  <*/.  ( 1985a). 

The  pertinent  equations  are  given  in  sections  II  and  III,  and  the  following 
constants  are  assumed: 


p  =  1025  Kg/m3 
pa=  1.2  Kg/m3 

a  =  2x  10~4  C_1 
Cp  =  3890  J  Kg'1  C_1 
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g  =  9.8  m/s2 
0  =  7.29 x  )0“5  s-1. 


For  each  of  the  above  defined  depths  and  for  each  latitude,  the  tunning 
constants  Cj  and  C2  will  be  adjusted,  so  that  the  summation  of  the  squares 

of  the  differences  between  the  diagnostic  mixed  layer  depths  and  the  observed 
mixed  layer  depths  is  a  minimum,  and  Cji  0  and  0  i  C2  i  2.  Following  this 

tuning  process,  a  standard  deviation  will  be  then  computed,  according  to  the 
definition: 


s2  =  1  (Hrhj)2/(N-1),  (4.4) 

»*  t 

where  the  Hj's  are  the  observed  mixed  layer  depths,  the  hj's  are  the 
diagnostic  mixed  layer  depths  and  N  is  the  number  of  longitude  points.  The 
standard  deviations  will  serve  as  a  quantitative  verification  of  the  models, 
the  model  assumptions,  and  the  accuracy  of  the  boundary  conditions. 

As  discussed  earlier,  the  rotation-stress  variable  is  required 
theoretically  to  be  greater  than  -1/2,  the  value  corresponding  to  a  situation 
of  zero  dissipation.  For  that  reason,  0  will  be  set  equal  to  -1/2  whenever 
its  calculated  value  is  less  than  -1/2.  During  the  month  of  September,  this  is 
likely  to  occur  in  the  western  Atlantic,  where  the  zonal  surface  wind  stress 
is  large  and  the  surface  heat  flux  is  small. 

The  selection  of  the  factor  R  and  the  penetration  depth  X  is  an 
imprecise  process  because  of  the  paucity  of  data  for  the  optical  properties  of 
the  equatorial  and  tropical  oceans.  Although  the  choice  here  of  values  for  R 
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and  X  was  based  on  a  preliminary  tunning  of  the  model  at  a  specific  site,  it  is 
still  somewhat  arbitrary.  Without  a  more  complete  data  set,  the  inclusion  of 
these  parameters  as  spacially-varying  tuning  parameters  does  not  seem 
justifiable  at  this  time.  Such  a  procedure  would  tend  to  mask  the  effects  of 
other  important  physical  mechanisms,  such  as  horizontal  and  vertical 
advection.  An  alternative  might  be  to  define  two  different  scales  of 
penetration  for  the  long  and  short-wave  radiation,  and  then  to  use  the 
experimental  results  of  Paulson  and  Simpson  (1977).  However,  this  would 
also  complicate  the  problem.  Furthermore,  Paulson  and  Simpson's  data  are 
valid  only  for  a  specific  spatial  and  temporal  situation  in  the  North  Pacific.  In 
the  absence  of  detailed  measurements  of  the  ocean's  optical  properties, 
there  is  no  simple  solution  to  this  problem,  and  the  results  of  the  present 
work  will  be  somewhat  affected. 


41 


C.  RESULTS 

Appendix  1  shows  the  zonal  variation  of  the  nondimensional  depth  H*  for 
latitudes  0,  10,  20  and  30  degrees  north.  Only  two  particular  cases,  the 
equator  and  10  degrees  north,  will  be  presented  and  discussed  in  the  text. 

1.  Equator 

Fig.  7  shows  the  equatorial  boundary  conditions  :  the  surface  heat  flux 
Q0,  the  net  surface  solar  irradiance  Qs,  the  total  surface  wind  stress  T,  and 
the  zonal  wind  stress  Tx.  Fig.  8a  shows  the  zonal  variation  of  the 
nondimensional  depth  H*.  computed  using  equation  (4.3).  It  is  apparent  from 
this  result  that  only  a  limited  zone,  between  about  31  degrees  west  and  4 
degrees  west  seems  to  be  in  a  steady-state  balance.  For  some  of  the 
remaining  regions,  the  presence  of  an  upward  total  surface  heat  flux 
precludes  the  possibility  of  an  equilibrium,  causing  the  top  of  the  thermocline 
to  continually  entrain.  We  have  previously  suggested  that  the  turbulent 
processes  near  the  eastern  boundary  are  likely  to  be  dominated  by  equatorial 
upwelling,  which  can  reduce  significantly  the  vertical  extent  of  mixing.  The 
values  from  Fig.  7a  seem  to  confirme  this  hypothesis,  since  the  values  of  H* 

are  less  than  1  near  the  coast.  For  the  western  region,  other  mechanisms 
might  partially  explain  the  departure  of  the  model's  results  from  the 
observations.  The  first  of  those  is  related  to  the  strong  wind-driven  current 
near  the  coast  of  Brazil,  the  "Guiana  current".  This  flow  and  the  associated 
meridional  variations  in  the  surface  heat  flux  might  cause  the  local  dynamics 
to  be  dominated  by  horizontal  advection.  A  second  possibly  contributing 
factor  is  the  influence  of  the  Amazon  River  runoff.  According  to  the  data  on 
the  mean  surface  salinities  in  Robinson’s  atlas, this  runoff  extends  to  more 
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Fig.  7(cont.).  Mean  surface  boundary  conditions  at  the  equator,  during 
September,  (c)  Total  surface  wind  stress,  (d)  Zonal  surface  wind  stress. 
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than  300  nautical  miles  eastward  from  the  coast  and  It  affects  the  upper  few 
tens  of  meters  of  the  water  column.  As  examined  by  Garvlne  (1984),  an 
advected  shallow  plume  of  buoyant  water  can  decrease  the  mixed  layer 
depth.  Also,  since  the  water  advected  from  the  Amazon  Is  more  turbid  than 
the  average  sea  water,  the  value  chosen  for  the  penetration  scale  X  is 
probably  too  large  for  the  region,  causing  the  model  to  overestimate  the 
mixed  layer  depth.  Fig.  7a  shows  precisely  this  effect  :  approaching  the 
western  boundary,  H*  falls  markedly  to  about  0.4. 

For  the  equatorial  region  between  31  degrees  west  and  4  degrees  west, 
the  diagnostic  values  were  calculated  for  the  mixed  layer  depths  scales  h0, 
h|  and  h2  .  The  results  were  compared  with  the  observed  mixed  layer  depths. 

Fig.  8  shows  the  zonal  profiles  obtained  for  each  of  the  defined  scales.  The 
vertical  bars  represent  the  effect  of  varying  Q0  by  ±10  W/m2,  which 

indicates  how  sensitive  the  model  is  to  moderate  errors  or  variations  in  the 
heat  flux  data.  The  comparaison  of  Figures  7a  and  8  also  shows  that  h2  is 

significantly  more  sensitive  to  variations  in  Q0  than  are  the  other  depth 
scales,  especially  when  Q0  is  small.  In  determining  a  best  fit  for  h0  and  h] , 
the  corresponding  standard  deviations  were  s0=  5.4  meters  and  sj=  6.2 
meters.  These  optimized  values  were  obtained  with  Cj=  1.3  and  Cj=  1.0,  for 
h0andhj,  respectively. 

Fig.  9  illustrates  the  tunning  process  of  Cj  and  C2,  in  computing  the 
zonal  variation  of  the  depth  scale  h2.  Displayed  here  are  isopleths,  in  m2,  of 

the  sums  of  the  squares  of  the  differences  between  the  diagnostic  and  the 
observed  depths,  as  a  function  of  C|  and  C2  values.  For  each  grid  point  of 

Ci  and  C2  values,  a  zonal  profile  of  the  diagnostic  depth  f»2  was  calculated, 
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Fig.  8.  Model  results  for  the  equator,  (a)  Zonal  variation  of  the 
nondimensional  mixed  layer  depth  H».  (b)  Zonal  variation  of  the  Kraus  and 

Turner  diagnostic  mixed  layer  depth  h0.  The  circles  represent  observed 
vali  es,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±10  W/m2.  The 

shaded  areas  in  (a)  correspond  to  the  zones  where  a  one-dimensional, 
steady-state  balance  is  assumed  impossible. 
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Fig.  8(cont.).  Model  results  for  the  equator,  (c)  Zonal  variation  of  the 
Kraus  and  Turner  diagnostic  mixed  layer  depth  h).  (d)  Zonal  variation  of  the 

Garwood  diagnostic  mixed  layer  depth  h2.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±10  W/m2. 
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as  well  as  the  zonal  summation  of  the  squares  of  the  differences  between 
those  depths  and  the  observed  mixed  layer  depths.  The  arrows  in  Fig. 9 

indicate  the  optimal  values  of  the  constants,  corresponding  to  a  minimum 
standard  deviation  s?  *  3.3  meters. 

There  are  some  interesting  features  to  be  noted  in  the  results*. 

(i)  There  is  no  quantitative  improvement  in  the  quality  of  the  results 
for  the  ICraus  and  Turner  model  (h0)  when  the  effect  of  the  penetrating 
radiation  is  included  (h^).  This  is  explained  by  the  fact  that  the  surface  heat 
flux  is  relatively  large  in  that  region,  and  the  vertical  scale  X  of  penetration 
is  much  less  than  the  computed  Obukhov  length  scale. 

(11)  The  zonal  variation  of  the  mixed  layer  depth  near  the  western  limit 
is  better  represented  when  the  rotation  stress  mechanism  is  included  (h?). 

In  that  region,  the  zonal  wind  stress  has  relatively  large  negative  values, 
which  causes  0  to  also  have  relatively  Urge  negative  values. 

(iii)  The  model  seems  much  more  sensitive  to  the  value  of  Cj  than  to 
the  value  of  C^.  Fig.  9  shows  a  relatively  large  zone,  between  C^= 1 . 2  and 
C2=l  .6  where  the  standard  deviation  remains  practically  constant.  This  fact, 
associated  with  the  similarity  of  the  results  between  h0  and  hj,  might  indicate 
that  the  effect  of  the  rotation  stress  mechanism  is  small, for  the  region  taken 
as  a  whole. 

In  general,  the  revised  model  shows  a  remarkably  good  agreement 
with  the  observations.  However,  and  because  the  Kraus  and  Turner  h0's  are 

already,  per  se,  in  good  agreement  with  the  observed  values,  the  results 
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are  not  conclusive  concerning  the  Importance  of  the  penetrating  radiation  and 
the  rotation  stress  mechanism.  The  eastern  part  of  the  equator, where  the 
observed  zonal  wind  stress  Is  positive,  would  constitute  a  suitable  domain  to 
further  test  the  significance  of  the  rotation  stress  mechanism.  There,  the 
diagnostic  depths  h?  would  be  larger  than  the  Kraus  and  Turner  depths  h0  and 

hf.  However,  as  we  have  seen,  other  mechanisms  appear  to  affect  the 
dynamics  of  mixing  In  the  region,  and  will  need  to  be  Included  in  such  a  test. 


Fig.  9.  Variation  of  the  sum  of  the  squares  of  the  differences,  in  m2  .between 
the  observed  equatorial  mixed  layer  depths  H  and  the  diagnostic  mixed  layer 
depths  hj  with  the  tunning  parameters  C|  and  C2. 
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2.  10  degrees  North 


Fig.  10  shows  the  mean  boundary  conditions  for  10  degress  north 

during  the  month  of  September.  In  Fig.  11a,  the  zonal  variation  of  the 
nondimensional  depth  H*  is  displayed.  This  figure  shows  that  the  zone  which 

seems  to  be  in  a  steady-state  balance  (between  about  38  degrees  west  and  17 
degrees  west)  is  considerably  smaller  than  the  corresponding  zone  for  the 
equator.  Although  the  surface  heat  flux  is  always  positive  for  the  entire 
domain,  the  values  are  usually  smaller  than  at  the  equator.  When  the  effect 
of  the  penetration  of  radiation  is  included,  the  result  is  to  further  reduce  the 
effective  heat  flux  at  the  surface,  as  defined  by  equation  (3.4).  This  effect 
reduces  buoyant  damping  and  precludes  the  possibility  of  an  equilibrium  for  a 
larger  area.  Although  unsteadiness  may  partially  account  for  these  results, 
advectlon  seems  to  be  a  more  plausible  explanation  for  the  regions  where  the 
effective  heat  flux  is  positive,  since  H*  is  always  less  than  1 .  A  meridional 

gradient  is  present  in  the  surface  heat  flux  because  of  the  strong  wind-driven 
current  in  the  western  part  of  the  tropical  Atlantic,  and  this  effect  might 

cause  buoyant  water  to  be  advected  northward,  reducing  the  vertical  extent 
of  mixing.  The  variation  of  H*  near  the  eastern  boundary  is  more  difficult  to 

explain.  From  about  18  degrees  west,  where  it  has  a  value  of  1 ,  H*  continues 

to  increase  eastward,  which  indicates  that  the  observed  depth  becomes 
larger  than  the  diagnostic  depth  h£.  The  analysis  of  the  temperature 

distribution  in  the  atlas  of  Robinson  et  at.  (1979)  shows  that  this  is  a  region 
of  particularly  strong  coastal  upwelling.  However,  positive  vertical 
advection  can  only  reduce  the  depth  of  the  turbulent  boundary  layer,  as 
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Fig.  10.  Mean  surface  boundary  conditions  at  10  degrees  north,  during 
September,  (a)  Net  downward  surface  heat  flux,  (b)  Net  downward  surface 
solar  irradiance. 


Fig.  10(cont).  Mean  surface  boundary  conditions  at  10  degrees  north,  during 
September,  (c)  Total  surface  wind  stress,  (d)  Zonal  surface  wind  stress. 
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stated  by  Garwood  et  al.  (1985b).  Thus,  this  mechanism  cannot  be 
responsible  for  the  difference. 

The  diagnostic  depths  h0,  h^  and  f>2  were  calculated  for  the  region 

between  38  degrees  west  and  17  degrees  west.  Fig.  11  shows  their  zonal 
variations.  The  vertical  bars  represent  the  effect  of  varying  Qo  by  ±5  W/rv  , 
from  which  it  can  be  concluded  that  the  model  is  much  more  sensitive  here  to 
small  variations  in  the  heat  flux  than  it  was  for  the  equator,  especially  when 
the  penetration  of  radiation  effect  is  included.  This  is  explained  by  the 
already  discussed  fact  than  the  surface  heat  flux  at  10  degrees  north  is 
significantly  smaller  than  it  is  at  the  equator.  In  some  areas,  particularly  in 
the  western  half  of  the  domain,  the  sensitivity  is  so  large  that  the  resulting 
depths  vary  by  a  factor  of  three  when  those  effects  are  considered. 
Considering  that  the  relative  errors  in  the  surface  heat  flux  data,  as 
estimated  by  Hastenrath  and  Lamb  (1978),  are  less  than  10  W/m2,  only  a 
qualitative  interpretation  of  the  results  for  this  latitude  should  be  made. 

For  each  of  the  diagnostic  scales  h0,  h1  and  h2,  the  optimized  tuning 

constants  and  the  corresponding  standard  deviations  have  the  following 
values: 

h<) :  *2.6,  C2  =  0.0,  s  =  11.8m 

hi  :  Cj  -  0.8,  C2  =  0.0,  s  =  7. 1  m 
h2  :  C-j  =  0.8,  C2’*0.7,  s  =  5.9m. 

A  significant  Improvement  in  the  quantitative  quality  of  the  results  is 
apparent,  when  the  penetration  effect  is  included  (hi ).  Note  that  the  optimum 
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Fig.  1i.  Mode]  resuits  for  10  degreees  north,  (a)  Zonal  variation  of  the 
nondimensiona)  mixed  layer  depth  H*.  (b)  Zonal  variation  of  the  Kraus  and 
Turner  diagnostic  mixed  layer  depth  h0.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±5  W/mZ.  The 
shaded  areas  in  (a)  correspond  to  the  zones  where  a  one-dimensional, 
steady-state  balance  is  assumed  impossible. 
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Fig.  I  l(cont).  Model  results  for  10  degrees  north,  (c)  Zonal  variation  of  the 
Kraus  and  Turner  diagnostic  mixed  layer  depth  h?.  (d)  Zonal  variation  of  the 
Garwood  diagnostic  mixed  layer  depth  h2.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±5  W/m2. 
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tuned  value  of  Cl  was  reduced  from  2.6  to  a  more  plausible  0.8,  which 
approaches  the  value  found  for  the  equator.  In  general,  the  Inclusion  of  the 
penetration  effect  tended  to  homogenize  the  diagnostic  mixed  layer  depth 
field.  This  caused  a  better  agreement  between  the  results  and  the 
observations,  particularly  near  the  extremes  of  the  longitudinal  Interval.  As 
expected  from  the  analysis  of  the  boundary  conditions  for  the  zonal  wind 
stress,  the  western  part  of  the  domain  is  better  represented  when  the 
rotation  stress  mechanism  is  included.  Also,  a  small  improvement  is 
observed  near  22  degrees  west,  where  the  observed  mixed  layer  depth  has  a 

minimum.  This  is  explained  by  the  zonal  wind  stress  being  positive  in  that 
region.  When  the  eastern  boundary  is  approached,  Tx  becomes  more  and 

more  positive,  causing  the  diagnosed  value  for  h?  to  be  less  than  the 
observed  value. 

Fig.  11  illustrates  the  tuning  process  for  h2-  Again,  the  contour  lines, 

in  m2,  of  the  sums  of  the  squares  of  the  differences  between  the  diagnostic 
and  the  observed  depths,  are  displayed  as  a  function  of  and  C2.  The 

optimal  value  of  C1  =  0.8  has  not  changed  with  the  inclusion  of  the  rotation 
stress  mechanism.  This  tends  to  verify  the  merits  of  the  procedure,  in  the 
sense  that  no  other  physical  mechanisms,  absent  from  the  model's 
formulation,  are  being  compensated  for.  On  the  other  hand,  although  Fig.  11 
shows  that  the  model  is  still  more  sensitive  to  variations  in  Cj  than  to 
variations  in  C2,  the  optimal  value  of  C2  is  much  better  defined  here  than  at 

the  equator.  Also,  this  relative  sensitivity  seems  to  verify  the  importance  of 
the  rotation  effect. 
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Although  a  much  better  representation  of  the  mixed  layer  depth  field 
was  obtained  with  the  revised  model  at  10  degrees  north,  this  verification  is 
limited  by  the  uncertainty  in  the  boundary  conditions.  As  we  have  already 
noted,  the  model  is  so  sensitive  to  small  variations  in  the  heat  flux  data,  that 
no  detailed  interpretation  of  the  results  is  justifiable. 


Fig.  12.  Variation  of  the  sum  of  the  squares  of  the  differences, in  m^, between 
the  observed  mixed  layer  depths  H  and  the  diagnostic  mixed  layer  depths  h2 

with  the  tunning  parameters  Cj  and  C2  (10  degrees  north). 
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Starting  with  the  theoretical  frame  work  of  the  one-dimensional, 
steady-state  Kraus  and  Turner  ( 1967)  model  for  the  surface  ocean  boundary 
layer,  a  revised  theory  was  presented.  This  new  theory  includes  the  effect  of 
the  penetrating  radiation  below  the  surface  on  the  buoyancy  flux,  as  well  as 
the  rotation  stress  mechanism  examined  by  Garwood  ( 1985a, b).  In  the 
presence  of  a  downward  surface  heat  flux  and  a  surface  wind  stress,  the 
penetration  of  radiation  will  tend  to  reduce  the  buoyant  damping  of 
turbulence,  enhancing  vertical  mixing.  The  rotation  stress  mechanism, 
driven  by  the  interaction  between  the  meridional  component  of  the  planetary 
rotation  Oy  and  the  zonal  wind  stress  Tx,  predicts  the  equilibrium  mixed 

layer  depth  to  increase  when  tx  is  negative  (westward)  and  to  decrease  when 
Tx  is  positive  (eastward). 

To  test  its  response,  and  to  gain  insight  into  the  relative  importance  of  the 
physical  mechanisms  involved,  the  revised  model  was  applied  to  a  set  of 
boundary  conditions  in  the  tropical  Atlantic.  Except  for  the  surface  wind 
stress,  where  digitized  values  were  used,  the  other  boundary  and  oceanic 
conditions  (observed  mixed  layer  depth,  surface  heat  flux  and  surface 
downward  solar  irradiance)  were  interpolated  from  low-resolution 
climatological  atlases.  Since  the  model  is  particularly  sensitive  to  small 
variations  in  the  heat  flux  data,  the  quality  of  the  results  is  considered  to 
have  been  somewhat  affected  by  the  amount  of  uncertainty  in  the  boundary 
conditions.  Also,  in  the  absence  of  detailed  measurements  of  the  optical 
properties  of  the  tropical  oceans,  the  values  of  the  penetration  parameters  X. 


and  R  were  assumed  to  be  constant  over  the  whole  domain.  Given  these 
limitations,  no  detailed  interpretation  of  the  results  is  appropriate  here. 

In  general,  the  response  of  the  revised  model  compares  favorably  with 
the  observations.  The  improvement  of  the  results  when  the  two  new 
mechanisms  are  included  is  very  significant,  particularly  at  10  degrees 
north.  This  was  expected,  in  so  far  as  the  penetration  effect  is  concerned, 
since  at  10  degrees  north  the  mixed  layer  is  relatively  shallower  than  at  the 
equator,  and  the  downward  surface  heat  flux  is  generally  smaller.  For  this 

latitude,  an  important  result  which  seems  to  support  the  need  to  include  the 
rotation  stress  formulation  is  that  the  optimal  value  of  the  constant  C] 

remains  constant  when  that  process  is  included  in  the  tuning  of  the  model. 

In  a  recent  paper,  Garwood  at  af.  (1985b)  examined  the  zonal  dependence 
of  the  climatological  mixed  layer  depth  In  the  near-equatorial  Pacific  and 
concluded  that  the  rotation  stress  mechanism  Is  a  plausible  explanation  for 
the  deep  mixing  In  the  central  equatorial  Pacific.  Although  that  feature  is  not 
so  pronounced  in  the  tropical  and  near-equatorial  Atlantic,  the  results 
obtained  here  seem  to  support  the  same  conclusion,  especially  In  Its  western 
part.  On  the  other  hand,  the  tropical  Atlantic  is  a  relatively  smaller  ocean 
basin,  and  the  influence  of  coastal  physical  processes,  like  coastal 
upwelling,  are  more  likely  to  Influence  significantly  the  dynamics  of  the 
boundary  layer.  Indeed,  this  effect  is  apparent  from  the  results,  since  only 
the  central  part  of  the  domain  seems  to  be  in  a  one-dimensional  equilibrium 
state.  Also,  the  sensitivity  of  the  model  to  small  variations  in  the  surface 
heat  flux  shows  that  this  equilibrium  state  is,  as  it  is  in  the  tropical  Pacific, 
a  delicate  one. 
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The  choice  of  the  tropical  Atlantic  as  a  domain  for  application  of  the  model 
was  primarily  motivated  by  the  existence  of  a  suitable  set  of  surface 
boundary  condition  data.  Future  research  on  the  subject  should  extend  the 
domain  to  the  non-tropical  regions  and  should  Include  unsteadiness  and 
advection.  Although  the  model  derived  in  this  paper  is  a  simplistic  one- 
-dimensional  representation  of  a  complex  phenomenon,  and  the  quality  of  the 
results  were  probably  affected  by  the  lack  of  detail  of  the  available  data,  it 
is  concluded  here  that  the  physical  mechanisms  of  rotation  stress  and 
penetration  of  radiation  are  important  in  determining  a  steady-state 
equilibrium  depth  of  turbulent  mixing  for  the  tropical  Atlantic. 
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Fig.  13.  Zonal  variation  of  the  nondimensionai  mixed  layer  aeptn  rt*.  iop  . 

equator.  Bottom  :  10  degrees  north.  The  shaded  areas  correspond  to  the 
zones  where  a  one-dimensional,  steady-state  balance  is  assumed 
impossible. 
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Fig.  14.  Zonal  variation  of  the  nondimensional  mixed  layer  depth  H*.  Top  :  20 

degrees  north.  Bottom  :  30  degrees  north.  The  shaded  areas  correspond  to 
the  zones  where  a  one-dimensional,  steady-state  balance  is  assumed 
impossible. 
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APPENDIX  B  -  COMPUTER  PROGRAMS 


'  Program  depth  -  calculates  the  steady-state  mixed  layer  depth 

variation  for  a  given  parallel  of  latitude. 

'  Jo  squint  Filips  Gaspar  -  April  1985  -  Naval  Postgraduate  School. 

DEFINT  i-k,n 

DIM  1  ong(200 ),  Qo(  200),  Qs  ( 200 ),  ta  u(  200 ),  taux(  200),  depth  ( 200  ),h(  200) 

********************  Variables  ******************************* 

lat  -  latitude  in  degrees  ;  long  -  longitude  in  degrees 

'  Qo  -  surface  heat  flux  in  W/m2  ;  Qs  -  surface  solar  irradiance  in  w/m2 
'  Q  -  effective  heat  flux  in  W/m2 

'  tau  -  total  wind  stress  in  N/m2  ;  taux  -  zonal  wind  stress  in  N/m2 
depth  -  observed  mixed  layer  depth  in  meters 
h  -  diagnostic  mixed  layer  depth  in  meters 
'  L  -  Obukhov  length  scale  in  meters 
Fi  -  rotation  stress  parameter 
'  R,lb  -  radiation  penetration  parameters 
*  cl  ,c2  -  model  tunning  parameters. 

************************************************************** 

Inputdata: 

INPUT  *  Enter  latitude  *;  lat :  lat$=STR$(  lat) 
fl$=*sh’+1at$  :  OPEN  *in*,*l,fl$ 

f2$=”qs'+latS  :  OPEN  ‘in",#2, f2$ 
f3$="tau’+lat$  :  OPEN  *in*,#3,f3$ 
f4$='taux'+lat$  :  OPEN  *tn*,M,f4$ 
f5$="mld'+lat$  :  OPEN  ’in",*5,f5$ 

n=0 

WHILE  NOT  EOF(l) 
n=n+l 

INPUT«l,long(n),Qo(n) 

INPUT*2,lg,Qs(n) 

IIIPUTi3,lgttau(n) 

INPUT»4,lg,taux(n) 

WPUT*5,lg,depth(n) 

WEND 

CLOSE 
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Inputparameters: 

PRIHT  ‘Longitude  limits  in  files  are  longO  ),long(n) 

INPUT  ‘Enter  new  longitude  limits  longO.longl 
IF  longO<1ong(l )  OR  long  1  >long(n)  THEN  Inputparameters 
INPUT  ‘Enter  value  of  Cl  cl 
INPUT  ‘Enter  value  of  C2  c2 
10  INPUT  ‘Enter  value  of  R  *;  R  :  IF  R<0  OR  R>1  THEN  10 

20  INPUT  ‘Enter  value  of  lambda  lb  :  IF  lb<=0  THEN  20 

Functions: 

DEF  FNL(tau,Q)=  124078!*tau*  1 . 5/Q 

DEF  FNQ(Qo,Qs,L)=Qo-R*Qs*(2*lb/L-(1+2*lb/L)*EXP(-L/lb)) 

DEF  F NFi(taux,Q)=1 44. 6*taux*C0S(  1  at/57. 2958  )/Q 

Calculate: 

sum=0 

FOR  1=1  TO  n 

IF  long(i)<1ong0  OR  long(i)>long1  THEN  30 
IF  Qo(i)<=0  THEN  h(i)=0  :  c$=‘*‘ :  GOTO  display 
L=FNL(tau(i),Qo(i)) 

Q=FNQ(Qo(i),Qs(i),L):IF  Q<=0  THEN  h(i)=0:c$=‘*‘:GOTO  display 
L=FNL(tau(i),Q) 

Fi=FNFi(taux(i) , Q)  :  IF  FK-.5  THEN  FI— .5 
h(i)=c1*L/( l+c2*Fi) 
sum=sum+(depth(i)-h(i))/'2 

c$-‘  ‘ 
display: 

PRINT  long(i);c$,h(i) 

30  NEXT  i 

PRINT  'Sum  :  *;  sum 

Choose: 

INPUT  ‘  Save  results(  1  ),New  fil es ( 2 ) ,  new  parameters(3),  Quit(4)  o 
ON  o  GOTO  Saveresults,  Inputdata,  Inputparameters, quit 

Saveresults: 

INPUT  ‘Enter  name  of  output  file  ‘;  fo$ 

OPEN  *out‘,»5,foS 
FOR  i=  1  TOn 

WRITE*5,long(i)  ,h(i) 

NEXT  i 
CLOSE 
GOTO  Choose 


*  Program  Spline  -  fits  a  cubic  spline  over  a  set  of  (x,y)  points. 

'  Joaquin  Filipe  Gaspar  -  January  1985  -  Naval  Postgraduate  School 

DIM  x(200),y(200),a(200,4),s(200)  :  DEFINT  i,j,n,e 
DIM  a0(200),al (200),a2(200) , a3(200) 

Enterdata: 

INPUT  ‘Read  disk  file(d)  or  enter  data  from  keyboard(k)  ";  aS 
IF  a$=*k*  THEN  key  ELSE  IF  a$="d"  THEN  disk  ELSE  Enterdata 

disk: 

INPUT  "Enter  name  of  file  to  fit  f$ 

OPEN  "in",*1,f$ 

n=1 

WHILE  NOT  EOF(  1 ) 

INPUT*  1  ,x(n)  ,y(n) 
n=n+1 

WEND 

CLOSE* 1 

n=n-1  :  GOTO  Endc 
key: 

INPUT  "Enter  number  of  points  n 

FOR  i=1  TO  n  :  PRINT  i,  :  INPUT  x(i),y(i)  :  NEXT  i 

Endc: 

INPUT  "Enter  end  condition  (1,2  or  3)  ";  endcond 
nm2=n-2  :  nm1=n-1 

dx1=x(2)-x( 1 )  :  dy1=(y(2)-y( 1 ) )/dx1*6 

defmatrix: 

FOR  1-1  TO  nm2 
dx2=x(i+2)-x(i+1 ) 
dy2=(y(i+2)-y(i+1 ))/dx2*6 
a(1, 1 )=dx1  :  a(i,2)=2*(dx1+dx2) 
a(i,3)»dx2  :  a(i,4)-dy2-dy1 
dx1=dx2 :  dy1=dy2 
NEXT  i 

ON  endcond  GOTO  solve ,  end2 ,  end3 
end2: 

a(i,2)=a(i,2)+x(2)-x(f) 
a(nm2,2)=a(nm2,2)+x(n)-x(nm! ) 

60T0  solve 
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end3: 

dx1»x(2)-x( 1 )  :  dx2*x(3)-x(2) 
a( 1 ,2)=(dx1+dx2)*(dxl+2*dx2)/dx2 
a( 1 ,3)=(dx2*dx2-dxl*dxl )/dx2 
dxn2-x(nm1  )-x(nm2)  :  dxnl*x(n)-x(nm1 ) 
a(nm2, 1 )=(dxn2*dxn2-dxn1*dxn1  )/dxn2 
a  ( nm2 , 2 )= ( dxn 1 +dxn2 )* ( dxn 1 +2*  dxn2 )/dxn2 

solve: 

FOR  i=2  TO  nm2 

a(i,2)=a(i,2)-a(i,  1  )/a(i-1 ,2)*a(i-1  ,3) 
a(i,4)=a(i,4)-a(i,  1  )/a(i-l  ,2)*a(i-l  ,4) 

NEXT  i 

a(nm2,4)=a(nm2,4)/a(nm2,2) 

FOR  i=nm2-l  TO  1  STEP  -1 

a(i,4)=(a(i,4)-a(i,3)*a(i+l  ,4))/a(i,2) 

NEXT  t 

FOR  i*1  TO  nm2 
s(i+1  )=a(i,4) 

NEXT  1 

ON  endcond  GOTO  endl 1 ,end22,end33 
endl  1: 

s(  1  )=0  :  s(n)=0:  60T0  coeff 
end22: 

s(  1  )=s(2)  :  s(n)=s(nm1 )  :  GOTO  coeff 
end33: 

s(  1  )=((dxHdx2)*s(2)-dx1*s(3))/dx2 
s(n)=(  (dxn2+dxn1  )*s(nm1  )-dxn1*s(nm2))/dxn2 

coeff: 

FOR  i=l  TO  nml 
dxl»x(i+1  )-x(i) 
a3(i)=(s(i+t )-s(i))/(6*dx1 ) 
a2(i)=s(i)/2 

al (i)=(y(i+l  )-y(i))/dx1-(2*dx1*s(i)+dxl*$(i+1  ))/6 
aO(i)=y(i) 

NEXT  i 

calculate: 

INPUT  "Enter  interval  betwen  points  dx 
PRINT  "x-limits  on  file  were  :  ",x(  1  ),x(n) 

INPUT  "Enter  new  x-limits  ";  xl  ,x2 
INPUT  "Enter  name  of  output  file  *;  f$ 

OPEN  "out*,#2,f$ 
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FOR  1=1  TO  nml 
xx-xl-x(i) 

WHILE  xl<x(i+l ) 

yy=aO(  f )+  a  1  ( i )  *xx+  aZC  i )  *xx*xx+  a3  ( i )  *xx*  xx*  xx 
WRITE*2,xl  ,yy 
xx=xx+dx :  x1=xl+dx 

WEND 
NEXT  i 

xx=x1-x(nml ) 

WHILE  xK=x2 

yy=a0(nm1  )+a1(nm1  )*xx+a2(0*xx*xx+a3(0*xx*xx*xx 

WRITE#2,x1,yy 

xx*xx+dx :  xl=x1+dx 

WEND 
CLOSE«2 
GOTO  Enter  data 


'  Program  Graph  -  plots  an  internally  defined  function  or  a  file  from 

disk. 

Joaquin  FUipo  Gaspar  -  Jan  1985  -  Naval  Postgraduate  School. 

DEFINT  a-z  :  DEFSNG  x,y,p  :  DIM  xf(700) ,y(700) ,yf(700) 
pi=3. 141593 :  YNF=3E+38 

menudef: 

MENU  1,0,1, ‘Function* 

MENU  1,1,1,  ‘Enter  new  function* 

MENU  1 ,2, 2, 'Use  internal  function’ 

MENU  1 ,3, 1 , ’Input  new  file  from  disk* 

MENU  1, 4,0, 'Use  present  file* 

MENU  2,0, 1 , ‘Resolution* 

MENU  2, 1,1, ’Very  High* 

MENU  2,2,1  ,*High‘ 

MENU  2,3,2, 'Medium* 

MENU  2, 4,1, -low- 
menu  3,0,1, ‘Scaling* 

MENU  3, 1,2, ’Automatic* 

MENU  3,2, 1 , ‘Enter  scaling  factors* 

MENU  4,0,1, 'X-Limits* 

MENU  4,1,1,  'Enter  x-limits* 

MENU  5,0,1, 'Grid* 

MENU  5, 1 ,2,*Yes* 

MENU  5,2,1,-No* 

MENU  6,0, 1 , "Title* 

MENU  6, 1,1, 'Enter  title* 

MENU  7,0,1,  "Run* 

MENU  7, 1,0, ‘Start  graph* 

MENU  7,2,0, ’New  graph* 

MENU  7,3,1, ’Stop-List* 

MENU  7,4,1, *Quit* 

MENU  ON 

re$=4  :  autoscale=1  :  grid=1 
Menuloop: 

IF  xrange  THEN  MENU  7,1,1 
ON  MENU  GOSUB  Menucheck 
GOTO  Menuloop 

Menucheck: 

ON  MENU(0)GOSUB  Function, Resolution, Scaling, Limits, Grid, Title, Start 

RETURN 


Function: 

ON  MENUO)  GOTO  Newfunction,01dfunction,Newfile,01dfile 

Newfunctlon: 

MENU  RESET  :  STOP 

Oldfunction: 

foption=0 

MENU  1,2,2:  IF  fflag  THEN  MENU  1,4,1 
RETURN 

No  wf  j  | 

INPUT  'Enter  name  of  file  *;  f$ 

OPEN  'in',«1,f$ 

i-0  :  yfmin-YNF  :  yfmax»-YNF 

WHILE  NOT  EOF(1) 

i=i+t  :  INPUT*  1,xf(1),yf(i) 

IF  yf(i)>yfmax  THEN  yfmax=yf(i)  :  imax=i 
IF  yf(f)<yfmin  THEN  yfm1n=yf(1)  :  1m1n=1 
WENO 
CLOSE* 1 

xfmfn=xf(  1 ) :  xfmax=xf(i) :  nx=1 
fflag*  1  :  foption=1  :  MENU  1,2,1:  MENU  1,4,2 
RETURN 
Oldfile: 

f option*  1  :  MENU  1,2,1  :  MENU  1,4,2 
RETURN 

Resolution: 

res=2~(MENU(1)-l) 

FOR  i*  1  TO  4  :  MENU  2, i,  1  :  NEXT  i 
MENU  2 ,  MENUC 1 ) , 2  :  RETURN 

Scaling: 

IF  MENU(  1  )=1  Tl'EN  autoscale*  1  :  MENU  3, 1 ,2  :  MENU  3,2, 1  '  RETURN 
autoscale-0  :  MENU  3, 1,1  :  MENU  3,2,2 
INPUT  'Enter  x-axis  step  ';  xstep 
INPUT  'Enter  minimum  y  ymln 
INPUT  'Enter  maximum  y  *;  ymax 
INPUT  'Enter  y-axis  step  ystep 
IF  (option  THEN  yrange=ymax-ymin 
RETURN 

Limits: 

ON  foption+1  GOTO  limitl  ,limit2 
limit  J: 

INPUT  'Enter  minimum  x  ';xmin  :  INPUT  'Enter  maximum  x  ';  xmax 
xrange*xmax-xmin 

IF  xrange*0  THEN  limit!  ELSE  RETURN 
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Iimit2: 

PRINT  'In  present  file,  minim  and  max  x’s  are  :  xfmin.xfmax 
INPUT  ’Enter  minimum  x  xmin  :  IF  xmin<xfmin  THEN  limit2 
INPUT  'Enter  maximum  x  xmax :  IF  xmax>xf  max  THEN  11mlt2 
xrange-xmax-xmin 

RETURN 


Grid: 

IF  MENU(  1  )=1  THEN  grid=1  :  MENU  5, 1 ,2  :  MENU  5,2, 1  :  RETURN 

grid=0  :  MENU  5,2,2  :  MENU  5,1,1  :  RETURN 


Title: 

t$="  :  INPUT  ‘Enter  title  *;t$  :  RETURN 
Start: 

IF  HENU(  1 )- 1  THEN  Graph  :  RETURN 
IF  ftENU(  1  )=3  THEN  MENU  RESET  :  STOP 
IF  MENU(  I  )-4  THEN  SAVE  'Graph' :  SYSTEM 
FOR  1-1  TO  6  :  MENU  1,0, 1  :  NEXT  i :  MENU  7,2,0 
CLS  :  CALL  TEXTSIZE02)  :  RETURN 


Graph: 

FOR  1-1  TO  6  :  MENU  1,0,0  :  NEXT  i 

CLS 

ON  f  option*  1  GOTO  graphl  ,graph2 

graph  1: 
i»0 

IF  autoscale=1  THEN  ymax=-YNF  :  ymin=YNF 
FOR  graphx=4  TO  488  STEP  res 
i=i+l 

x=xm1n+xrange*(graphx-4)/484 

GOSUB  Fx 

IF  autoscale=0  THEN  30 

IF  y(1)>ymax  AND  ABS(yO)oYNF)  THEN  ymax=y(i) 
IF  y(i)<ymin  AND  ABS(y(i)<>YNF )  THEN  ymw»y(i) 
30  NEXT  graphx 

yrange=ymax-ymin :  ni=i-l 
GOTO  setitle 

graph2: 

i=l 

WHILE  xf(i)<xm1n 
i«i+l 

WEND 

ii=i  :  xmin=xf(i) 
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'v» 


I 


a 


<=nx 

WHILE  xf(i)>xmax 
i=i- 1 

WEND 

ifn-i :  xmax-xf(i)  :  xrange-xmax-xmin 

IF  autoscale  THEN  ymax=yfmax  :  ymin=yfmin  :  yrange=ymax-ymin 
it  1 8* 

CALL  HOVETO(  100,20)  :  PRINT  t$; 

xlabel: 

CALL  TEXTSIZE(9) 

LINE(4,4)-(4,284)  :  LINE(4,284)-(488,284) 

IF  autoscale  THEN  x$tep=xrange/4 
FOR  x=xmin+xstep  TO  xmax  STEP  xstep 
graphx=-484*(xmin-x)/xrange+4 
PSET(graphx,283) 

IF  grid-0  THEN  40 
FOR  graphy=280  TO  4  STEP  -4 
PSET(graphx.graphy) 

NEXT  graphy 

40  CALL  MOVETO(graphx-12,295) 

IF  xoxmax  THEN  PRINT  x; 

NEXT  x 

ylabel: 

IF  autoscale  THEN  ystep=yrange/4 
FOR  y=ymin  TO  ymax  STEP  ystep 

graphy=280*(ymin-y)/yrange+284 
PSET( 5, graphy) 

IF  grid=0  THEN  50 
FOR  graphx=9  TO  486  STEP  4 
PSETtgraphx,  graphy) 

NEXT  graphx 

50  CALL  MOVETO(6,graphy-2) 

PRINT  y; 

NEXT  y 

Plot: 

ON  f option*  1  GOTO  plotl  ,plot2 
plot  1 : 

FOR  i=1  TO  ni :  IF  (y(i)>ymax  OR  y(i)<ymin)  THEN  NEXT  i 
gy=280*(ymin-y(i)  )/yrange+284 
gx-4+( i-1  )*re$  :  i-0 
FOR  graphx=4  TO  488  STEP  res 
i=i+l 
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IF  (y(0>ymax  OR  y(i)<ymin)  THEN  flagy=0  :  GOTO  5 
flagy*flagy+1  :  graphy=280*(ymirt-y(i))/yrange+284 
IF  f lagy>  1  THEN  LINE(gx,gy)-(graphx,graphy) 
gx=graphx :  gy=graphy 

NEXT  graphx 
MENU  7,2,1 

RETURN 

p1ot2: 

gx*-484*  ( xmin-xf  ( il )  )/xr  ange+4 
gy=280*(ymin-yf(ii))/yrange+284 
FOR  1=i1  TO  ifn 

graphx=-484*(xm1n-xf(i))/xrange+4 
graphy=280*(ymin-yf(i))/yrange+284 
LINE(gx,gy)-(graphx,graphy) 
gx=graphx :  gy=graphy 

NEXT  i 

MENU  7,2, 1 

RETURN 


y(i)=EXP(-x*x) 

RETURN 
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